clear all

set more off
set matsize 10000

global base "X:\Dropbox\educacion superior 2021\submission\EDCC\For publication\Replication package\Data"
global output "X:\Dropbox\educacion superior 2021\submission\EDCC\For publication\Replication package\Output"

cd "$base"
use "output/panel_final"

*===============================================================================
* Robustness - different definitions of skills
*===============================================================================

* Skills (+ respect in r1)

sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> respect_r1) ///
(skills_c_r1 respect_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 respect_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 order_1 msib_r3 fsib_r3 migrated_r1 indigena_r1 skills_cr_r3 skills_nc2_r3 -> enroll_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model

{
cap mat drop Rob_Skills1
mat Rob_Skills1 = J(8,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc2_r3 {
		
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat Rob_Skills1[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup:`s'.sex#c.`var']/_se[enroll_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat Rob_Skills1[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

* Skills (only validated measures)

sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc1_r2) ///
(skills_cr_r2 skills_nc1_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc1_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc1_r3) ///
(terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 order_1 msib_r3 fsib_r3 migrated_r1 indigena_r1 skills_cr_r3 skills_nc1_r3 -> enroll_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc1_r2 e.skills_cr_r3*e.skills_nc1_r3) method(mlmv)

beep

** storing coefficients of enrollment model

{
cap mat drop Rob_Skills2
mat Rob_Skills2 = J(8,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc1_r3 {
		
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat Rob_Skills2[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup:`s'.sex#c.`var']/_se[enroll_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat Rob_Skills2[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

*** Outputting results

{
qui putexcel set "$output/AppTable-RobustSkills.xlsx", replace

putexcel C4 = matrix(Rob_Skills1), hcenter
putexcel G4 = matrix(Rob_Skills2), hcenter

qui putexcel 	C2 = "A" 		///
				G2 = "B" 		///
				C3 = "Male"		///
				E3 = "Female"	///
				G3 = "Male"		///
				I3 = "Female"	///
				, hcenter

local v = 0
foreach var in terwi2_r3 terwi3_r3 cognitive noncognitive {
	qui putexcel B`=4+`v'*2' = "`var'"
	local ++v	
}
}

*===============================================================================
* Robustness - different definitions of enrollment
*===============================================================================

* Skills (still enrolled in sup)

sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 order_1 msib_r3 fsib_r3 migrated_r1 indigena_r1 skills_cr_r3 skills_nc2_r3 -> still_sup) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model

{
cap mat drop Rob_Enroll1
mat Rob_Enroll1 = J(8,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc2_r3 {
		
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat Rob_Enroll1[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[still_sup:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[still_sup:`s'.sex#c.`var']/_se[still_sup:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat Rob_Enroll1[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

* Skills (+ conservador: sabatico podria enroll)

sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 order_1 msib_r3 fsib_r3 migrated_r1 indigena_r1 skills_cr_r3 skills_nc2_r3 -> enroll_sup1) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model

{
cap mat drop Rob_Enroll2
mat Rob_Enroll2 = J(8,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc2_r3 {
		
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat Rob_Enroll2[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup1:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup1:`s'.sex#c.`var']/_se[enroll_sup1:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat Rob_Enroll2[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

* Skills (- conservador: 0 a todos los que no estan enrolled yet)

sem (migrated_r1 indigena_r1 careed_sec careed_ter agechild_r1 zhfa_r1 order_1 msib_r1 fsib_r1 rural_r1 terwi2_r1 terwi3_r1 -> skills_c_r1) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_cr_r2) ///
(skills_c_r1 parinv_r1 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r2 zhfa_r2 order_1 msib_r2 fsib_r2 rural_r2 terwi2_r2 terwi3_r2 -> skills_nc2_r2) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_cr_r3) ///
(skills_cr_r2 skills_nc2_r2 parinv_r2 migrated_r1 indigena_r1 zhfa_r1 careed_sec careed_ter agechild_r3 zhfa_r3 order_1 msib_r3 fsib_r3 rural_r3 terwi2_r3 terwi3_r3 -> skills_nc2_r3) ///
(terwi2_r3 terwi3_r3 careed_sec careed_ter rural_r3 zhfa_r3 order_1 msib_r3 fsib_r3 migrated_r1 indigena_r1 skills_cr_r3 skills_nc2_r3 -> enroll_sup2) ///
if panel1234==1, group(sex) vce(cl clustid) cov(e.skills_cr_r2*e.skills_nc2_r2 e.skills_cr_r3*e.skills_nc2_r3) method(mlmv)

beep

** storing coefficients of enrollment model

{
cap mat drop Rob_Enroll3
mat Rob_Enroll3 = J(8,4,.)

local i=1
foreach var in terwi2_r3 terwi3_r3 skills_cr_r3 skills_nc2_r3 {
		
	forvalues s=0/1 { // male then female
		
		local c = 1
		foreach coef in b se { // coef then se
			mat Rob_Enroll3[`=(`i'-1)*2+`c'',`=`s'*2+1'] = `=round(_`coef'[enroll_sup2:`s'.sex#c.`var'],0.01)'
			local c = 2
		}
		
		local pval = 2*(normal(-abs(_b[enroll_sup2:`s'.sex#c.`var']/_se[enroll_sup2:`s'.sex#c.`var'])))
		local stars `=cond(`pval'<0.01, 3, `=cond(`pval'<0.05, 2, `=cond(`pval'<0.1, 1, 0)')')'

		mat Rob_Enroll3[`=(`i'-1)*2+1',`=`s'*2+2'] = `stars'

	}
	local ++i
}
}

*** Outputting results

{
qui putexcel set "$output/AppTable-RobustEnroll.xlsx", replace

putexcel C4 = matrix(Rob_Enroll1), hcenter
putexcel G4 = matrix(Rob_Enroll2), hcenter
putexcel K4 = matrix(Rob_Enroll3), hcenter

qui putexcel 	C2 = "A" 		///
				G2 = "B" 		///
				K2 = "C" 		///
				C3 = "Male"		///
				E3 = "Female"	///
				G3 = "Male"		///
				I3 = "Female"	///
				K3 = "Male"		///
				M3 = "Female"	///
				, hcenter

local v = 0
foreach var in terwi2_r3 terwi3_r3 cognitive noncognitive {
	qui putexcel B`=4+`v'*2' = "`var'"
	local ++v	
}
}

*===============================================================================